Properties of Quantum Systems via Diagonalization of Transition 
Amplitudes II: Systematic Improvements of Short-time Propagation 
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In this paper, building on a previous analysis |l| of exact diagonalization of the space-discretized 
evolution operator for the study of properties of non-relativistic quantum systems, we present a 
substantial improvement to this method. We apply recently introduced effective action approach for 
obtaining short-time expansion of the propagator up to very high orders to calculate matrix elements 
of space-discretized evolution operator. This improves by many orders of magnitude previously used 
approximations for discretized matrix elements and allows us to numerically obtain large numbers 
of accurate energy eigenvalues and eigenstates using numerical diagonalization. We illustrate this 
approach on several one and two-dimensional models. The quality of numerically calculated higher 
order eigenstates is assessed by comparison with semiclassical cumulative density of states. 
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I. INTRODUCTION 

In first paper in this series [ll we analyzed in detail 
the earlier introduced method 

H H II m for studying 
properties of quantum systems based on the diagonaliza- 
tion of real-space discretized evolution operator, as well 
as the errors associated with the discretization process. 
This analysis provided us with a better understanding of 
this method that can be used to numerically calculate 
energy eigenvalues and eigenstates of few-body physical 
systems. We have shown that errors due to the finite 
discretization step A used for space discretization vanish 
exponentially with 1/A 2 for short time of propagation. 
This highly outperforms the usual (polynomial in A 2 ) 
behavior of errors in approaches with diagonalization of 
space-discretized Hamiltonians @, H B Q • 111 addition, 
derived analytic estimates for discretization errors pro- 
vide an easy way for taking into account such errors and 
eliminating them inpractical applications. For these rea- 
sons, the approach [2| now becomes the preferred method 
to study systems with few degrees of freedom, for which 
it can be efficiently and straightforwardly implemented. 
In our previous paper we have also analyzed estimates 
for errors due to the introduction of the space cutoff L, 
providing simple criterion for assessing and eliminating 
errors of this type. 

In the previous paper we have also demonstrated that 
time of propagation t, a parameter introduced by this 
method, is a source of a new type of error that comes 
about from using short time approximations. This prob- 
lem was not addressed at all in Ref. 0. It has recently 
been discussed [ltj HH, G3 and will be the main focus of 
the present paper. Errors associated with the time of evo- 
lution parameter t must be carefully taken into account 
and may substantially limit the precision of numerical 
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calculation in the diagonalization method. In this pa- 
per we address this problem by applying the recently 
introduced effective action approach [l]J Q, [HI, [lj| [lTj 
to systematically improve approximations for transitions 
amplitudes. This in turn leads to the many orders of 
magnitude decrease of errors in obtained energy eigen- 
values, as shown in this paper. We demonstrate on sev- 
eral lower-dimensional models how use of higher-order ef- 
fective actions significantly reduces numerical errors and 
systematically improves obtained energy eigenvalues and 
eigenstates. 

Together with the results from our previous paper, this 
paper completes the analysis of the method based on 
the diagonalization of transition amplitudes, providing us 
with necessary analytical knowledge to estimate errors of 
all types associated with this method and to numerically 
very accurately calculate large numbers of energy eigen- 
values and eigenstates. This invites various applications 
of the method to the study of few-body quantum systems, 
some of which are discussed throughout the paper. 

The text is organized as follows: in Section[TT]we briefly 
review the effective action approach and demonstrate 
how it can be used for numerical calculation of transition 
amplitudes. In Sections IIIII and IIVI we apply the exact 
diagonalization method [J, 0] improved by the use of ef- 
fective actions to the numerical study of several one and 
and two-dimensional models. In these sections we also 
show how the number of reliable energy eigenvalues can 
be estimated using comparison of numerically obtained 
results with semiclassical cumulative density of states for 
higher-lying eigenstates. Section fVl gives our concluding 
remarks and some relevant applications of this approach. 



II. EFFECTIVE ACTIONS 

To introduce the notation, we first briefly outline the 
diagonalization method Q, presented in more detail in 
Section 2 of our previous paper After discretizing 
the continuous space and replacing it with a grid defined 
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by a discretization step A, all the quantities are denned 
only on a discrete set of coordinates x n = nA, where 
n G Z is any integer number. For a physical system with 
Hamiltonian H, the evolution operator (in the imaginary 
time formalism) is defined as exp(-tH), where t is the 
time of evolution. Transition amplitudes are defined as 



A{x,y;t) = (x\e tH \y), 



(1) 



and give the discretized evolution operator matrix ele- 
ments 



Anm (0 



A(nA, mA; t) , 



(2) 



where d is the number of spatial dimensions. The 
eigenvectors of such a matrix correspond to the space- 
discrctized eigenfunctions of the original Hamiltonian, 
while the eigenvalues are related to the eigenvalues of 
the Hamiltonian and can be written as 



-tEk (A.L.t ) 



(3) 



where we emphasize the dependence of the numerically 
calculated eigenvalues on all discretization parameters. 
The number of obtained eigenvalues and eigenstates is 
equal to the linear size of matrix A, which has to be finite 
when we represent any physical system on the computer. 
Typically, we restrict the range of indices n, m to the 
finite interval — N < n,m < N, so that the number of 
points in the grid is S = (2N) d . Note that the range 
can be adjusted so that the size S is an odd number. In 
Eq. (|5J) we have also introduced the space cutoff L, which 
corresponds to the restriction on the range of grid-point 
indices n,m, and is given by L = NA. 

As we can see, the precise calculation of transition 
amplitudes is essential for practical applications of this 
method. In Ref. @ all calculations are based on the naive 
approximation for transition amplitudes 



AV>(x,y;t) 



1 



(27rf) d / 2 



(4) 



which yields energy eigenvalues correct only to order 
0(t), and is for this reason designated by A^\ If one 
uses the naive approximation for transition amplitudes, 
then times of propagation must be very short for errors 
to be small enough. Practically, even for short times of 
propagation, such errors are always much larger than the 
errors due to discretization, and therefore significantly 
limit the applicability of the method. In addition to this, 
the results obtained in our previous paper [l[ on exactly 
solvable models suggest that longer times of propaga- 
tion generally give smaller errors in the diagonalization 
approach. The trade-off between these effects and its im- 
plications on numerical results have been documented in 

To address thiSjin principle one can use Monte Carlo 
simulations [HI, Eil to calculate amplitudes A to high 
precision. Although this can effectively resolve the prob- 
lem in many cases, it is often numerically very expensive. 




FIG. 1: (Color online) Transition amplitude A^ p '(Q,Q;t) as a 
function of the time of propagation t, calculated analytically 
using different levels p of the effective action. The plot is for 



the quartic anharmonic potential V(x) 
with parameters M = u = 1, g = 10. 
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More importantly, resorting to the use of Monte Carlo 
practically limits further analytical approaches. We will 
instead use the recently introduced effective action ap- 
proach [3, 0, EH, Eil E3] that gives closed- form analytic 
expressions A^ (x, y; t) for transition amplitudes which 
converge much faster, 



A^(x,y;t)=A(x,y;t) + 0(tP +1 /t d / 2 



(5) 



where p is an integer number corresponding to the order 
of the effective action used, i.e. order of energy eigen- 
value errors t p . For a general many-body theory effective 
actions up to p = 10 have been derived, while for a spe- 
cific models much higher values can be obtained, e.g. for 
the anharmonic oscillator and other polynomial interac- 
tions, for which effective actions have been calculated up 
top = 144. So, if p is high enough, it is sufficient that the 
time of evolution is less than the radius of convergence 
of the above series (t < t c ~ 1) and errors in calculated 
values of transition amplitudes will be negligible. This is 
illustrated in Fig. [T] for the case of a quartic anharmonic 
oscillator. The use of high-order expansion in the time 
of propagation of amplitudes will allow us to use times of 
evolution up to r c , which are much longer than the typi- 
cal times one can use with the naive (p = 1) amplitudes. 
At the same time, the expansion up to very high orders 
substantially decreases the errors associated with t, and 
may practically eliminate them. 

The analytic expressions for higher-order approxima- 
tions for transition amplitudes are based on the notion 
of effective actions, which are introduced by casting the 
solution of the time dependent Schrodinger equation for 
the transition amplitude in the form 



A(x,y;t) = 



1 



(27rf) d / 2 



-tw( a - 



-y;t) 



(0) 
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where W(x, 5; t) is the effective potential, with the fol- 
lowing boundary behavior: 

timW(x,5;t) = V(x). (7) 

As shown previously, the effective potential W(x,5;t) is 
regular in the vicinity of t = 0, enabling us to represent it 
in the form of a power series in short time of propagation 
t. The coefficients in this series are functions of the po- 
tential and its derivatives. The truncation of the series 
for the effective potential up to order designated 
by wAp -1 ) (x, 5; t), gives the expansion of the transition 
amplitude accurate to t p+1 /t d / 2 , 

A^(x,y;t) = j^e-^-^^^. (8) 

The analytic expressions for higher-order effective actions 
therefore yield analytic approximations for amplitudes 
with the convergence behavior given by Eq. ([5]) . We em- 
phasize that although the structure of the effective action 
solution form ([6|) is motivated by the path integral for- 
malism, the expression for amplitudes obtained in the 
above approach contains no integrals and can be used 
straightforwardly as long as the time of propagation is 
below the radius of convergence of the short-time series 
expansion. 

For the exactly solvable case of a harmonic oscillator 
one finds that the radius of convergence is r c = tt/uj. The 
radius of convergence is simply the distance in the com- 
plex time plain from the origin to the nearest singularity 
of the propagator. For the harmonic oscillator the singu- 
larities are located at ±ikir/u), k G N. The consequence 
of these singularities is that the power series for the ef- 
fective potential W(x, S;t) converges only for t < t c . It 
is often difficult to analytically determine the radius of 
convergence of the short time expansion of the transition 
amplitude. However, numerically this is a very simple 
problem, since outside of the radius of convergence the 
calculated approximative amplitudes rapidly tend to in- 
finity (for levels p for which the effective potential is not 
bounded from below; see Ref. [2(|) or to zero with the in- 
crease of p. From Fig. [T] we easily estimate radius of con- 
vergence to be t c ~ 1 for a quartic anharmonic potential 
V(x) = \Mup-x 1 + ^x A , with parameters M = u = 1, 
g = 10. Such numerical determination of the radius of 
convergence for a given level p is always done before prac- 
tical use of the effective potential. Note that we are not 
interested in the precise value of r c , just in its rough 
estimate which will allow us to safely use times of prop- 
agation below r c . 

To conclude the section, let us stress that the effective 
action approach can be used only for sufficiently smooth 
potentials, i.e. those that have derivatives of the required 
order, corresponding to the level p of effective action, as 
discussed in Ref. [lj]. For potentials that do not fulfill 
this condition (e.g. stepwise potentials), the effective ac- 
tion approach cannot be directly used. However, one can 
replace the original potential with some of its smooth 



deformations, perform numerical calculations, and at the 
end take the limit of the deformation parameter in which 
the original potential is recovered. The numerical results 
obtained in such a way must be carefully cross-checked 
using other methods. 

III. NUMERICAL RESULTS FOR d = 1 MODELS 

In this section we apply the approach outlined above 
to several d = 1 models and demonstrate its substantial 
advantages for numerical studies of eigenstates of vari- 
ous physical systems. We numerically analyze all sources 
of errors present in this approach due to discretization 
parameters L and A, as well as the time of propagation 
parameter t. We present the obtained numerical results 
for energy eigenvalues and eigenstates. We also assess the 
quality of the obtained energy spectra through compari- 
son with the semiclassical approximation for the density 
of states, which should be accurate at least for the higher 
regions of the spectrum. 

The first model we study is the quartic anharmonic 
oscillator with potential 

V(x) = -Mtu 2 x 2 + -^-x 4 . (9) 

For this potential the effective actions have been previ- 
ously derived up to p = 144 [211 ] . and here we will use 
various levels p to illustrate the dependence of errors on 
the level p used in calculations. 

Fig. [U presents the analysis of various errors in the 
ground energy calculation for a particular choice of pa- 
rameters of the potential M = u> = 1, g = 48. The 
spectrum of the potential is calculated by the numerical 
diagonalization of the space-discretized transition ampli- 
tude matrix. The errors are estimated using the exact 
value of the ground energy calculated elsewhere [22j by 
a different technique to very high precision. The depen- 
dance of the error related to the introduction of the space 
cutoff L is illustrated in Fig. [2^, while Fig. [2Jd gives the 
dependence of ground energy errors on the time of prop- 
agation parameter t for various values of the discretiza- 
tion step A. On both graphs we see the results obtained 
with effective actions of different levels p. Fig. [2Jd clearly 
shows that the errors due to the time of propagation are 
proportional to t p , as expected when we use the effective 
action of the level p. The errors in eigenvalues are of the 
same order as errors in calculation of individual matrix 
elements, and for this reason we see the typical t p behav- 
ior of ground and higher energy eigenvalues. It is already 
now evident that the use of higher order effective actions 
increases the accuracy of numerically calculated energy 
eigenstates for many orders of magnitude. This is the 
most important contribution of this paper. 

The L-dependence of the error is analytically known 
0, 0] • The saturation of errors in Fig. [2^ for a given 
level p corresponds to a maximal precision that can be 
achieved with that p, i.e. denotes the value of L for 
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which errors introduced by other sources become larger 
than the error due to the finite value of the space cutoff. 
This can be easily seen if we combine the data from both 
graphs. For example, the level p = 9 effective action has 
the saturated value of the error of the order of 10 -14 . 
For t — 0.02 we find that the error due to the time of 
propagation is of the same order if one uses sufficiently 
fine discretization (A = 0.05). Therefore, the saturation 
of errors on the left-hand graph are caused by the errors 
due to the time of propagation. However, if one uses dis- 
cretization which is not sufficiently fine, the saturation 
of errors can be also caused by the discretization effects. 
Such effects can be also analytically estimated to be pro- 
portional to -2cxp(-27r 2 £/A 2 )cosh(7r 2 (fc + l)t/LA)/t, 
as we have shown in the previous paper [l[ . 

TableUgives low-lying energy eigenvalues of the anhar- 
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FIG. 2: (Color online) (a) Deviations from the ground energy 
\Eq P \A, L,t) — Eo xact \ as a function of the space cutoff L and 
(b) as a function of the time t. The ground energy is obtained 
using different levels p = 1,3,5,7,9,11,13 (top to bottom) 
of the effective action for the quartic anharmonic potential, 
with parameters M = uo = 1, g = 48, A = 0.05, t = 0.02 
on graph (a), and L — 4 on graph (b). The exact ground 
energy_E^ act = 0.95156847272950001114693 . . . is taken from 
Ref. [23]. Dashed lines on the graph (b) correspond to the 
known discretization error 



k 


Eh 




5 Ei, 





U.9515d84727295UU(J1114d8(8J 


5 x IO - '" 


6 x 10 _:ia 


1 


o.z\)zob ( ozl4o44Doyzzoyi(D I ) 


4 X 1U 


o w 1 a - 
z X 1U 


z 


D.0U000UOD ( 44DOZDUyyoy(OZ^ ) 


1 X 1U 


/I w 1 A - 

4 X 1U 


6 


a 707000 1 70 7no7n rni r r / a a o\ 
6161 ( z fUofUDU100o(445j 


X 1U 


O X 1U 


A 

4 


lo.45lZ /OOODUoOOOyoOOO(14c5y ) 


A X 1U 


z X 1U 





1 / .oi4iozoyyzoou / uyzoy i^dzud j 


X 1U 


z X 1U 


6 


21.7909563917965158973(8744) 


6 x 10"™ 


3 x lQ- 2i 


7 


26.286125156056810490(92289) 


2 x 10~ ia 


7 x lO^ 1 


8 


30.979882837938369575(08213) 


2 x 10- ia 


8 x lO^ 1 


9 


35.856438766665971146(24181) 


3 x 10" iy 


9 x 1Q- Z1 



TABLE I: Low-lying energy levels of the anharmonic quartic 
potential, obtained by diagonalization using level p — 13 effec- 
tive action. The parameters are M = u) = 1, g = 48, L — 5, 
A = 0.05, t — 0.01. For higher energy eigenvalues, absolute 
and relative errors AEk and SEk are estimated by comparison 
with the diagonalization results obtained from higher-order ef- 
fective actions, finer discretizations, larger space cutoffs, and 
lower values of the propagation time t. 



monic oscillator for a particular choice of the parameters 
of the potential and discretization parameters. In princi- 
ple, one can achieve arbitrary precision by the use of ap- 
propriately chosen discretization parameters. Of course, 
for arbitrary precision calculations one has to use one of 
the software packages able to support such calculations. 
For example, we have used Mathematica [25| in order to 
be able to achieve high-precision results presented on the 
above graphs. The important conclusion is that even for 
very moderate values of discretization parameters, the 
use of higher-order effective actions leads to very small 
errors, which may be practically implemented with min- 
imal computing resources. 

The analysis of errors such as the one presented in 
Fig.[2]is sufficient to estimate optimal values of discretiza- 
tion parameters. In general, for a desired numerical pre- 
cision of energy eigenvalues, the optimal values of param- 
eters are chosen so that all types of errors are approxi- 
mately the same. The overall error is always dominated 
by the largest of all errors, and therefore it is optimal to 
have all errors of the same order of magnitude. 

For specific calculations one can have additional con- 
straints. For example, if one is interested only in energy 
eigenvalues, then the optimal parameters are obtained by 
minimizing all errors and minimizing the ratio N — L/A, 
which corresponds to the size of the transition operator 
matrix S = 2N that needs to be numerically diagonal- 
ized. The minimization of N is performed in order to 
minimize computation time needed for the diagonaliza- 
tion, which roughly scales as iV 3 . 

On the other hand, if one is interested in details of en- 
ergy eigenf unctions, then it might be necessary to have 
a fixed small value for the discretization step A, which 
will allow all features of eigenstates to be visible. This is 
especially important for studies of higher energy eigen- 
functions which e.g. have many nodes, and in order to 
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TABLE II: Low-lying energy levels of the double-well poten- 
tial, obtained by diagonalization using level p = 18 effective 
action. The parameters used: M = — 1, u> = 1, g = 12, 
L = 16, A = 0.1, t — 0.05. The absolute and relative errors 
AEk and SEk are estimated by comparison with the diago- 
nalization results obtained from higher-order effective actions, 
finer discretizations, larger space cutoffs, and lower values of 
the propagation time t. 



density of states, defined formally as 



p{E) =Y j 8{E - E k 



(10) 



k=0 



for a system with a discrete spectrum. This relevant 
physical quantity can be directly calculated using nu- 
merically obtained spectra. On the other hand, it can be 
also analytically calculated using semiclassical approxi- 
mation. This approximation is valid at least in the high- 
energy region, and we can use it to assess the quality 
of our numerical results. In semiclassical approximation, 
the density of states is calculated as 



Psc{E) 



d d xd d p 

(2irh) d 



6{E - H(x,p)) . 



(11) 



replacing the discrete spectrum with a continuous distri- 
bution of energy defined by the classical Hamilton func- 
tion -ff (x, p). After integration over momenta, we obtain 
the well known result 123 



study them it is necessary to have sufficient spatial res- 
olution. In such case, the value of A is fixed and other 
parameters are chosen so as to minimize the errors to a 
desired value. For example, with the discretization step 
of the order A = 10 -3 we have been able to accurately 
calculate several hundreds energy eigenfunction of the 
quartic anharmonic oscillator. 

Table [TT] gives eigenvalues of the double- well potential, 
obtained from the quartic anharmonic potential §§§ by 
setting the mass M to some negative value. As can be 
seen, numerically obtained energy eigenvalues have the 
precision similar to the previous case of the quartic po- 
tential without symmetry breaking. The double well be- 
havior of the potential does not present any obstacle in 
its numerical treatment by this method. 

Another situation in which one might be interested 
to keep the ratio N = L/A, i.e. the size of the space- 
discretized evolution operator matrix as large as possible 
is when a large number of energy eigcnlcvcls is needed. 
The number of energy eigenvalues that can be calculated 
by the diagonalization is limited by the size of the ma- 
trix S — 2N. Usually the highest energy levels cannot 
be used due to the accumulation of numerical errors, and 
therefore one needs to have a matrix of sufficient size in 
order to study energy spectra. In such cases it is neces- 
sary to use highly optimized libraries for numeric diago- 
nalization. We have implemented the effective actions as 
a C programming language code [2JJ and used LAPACK 
[2(| library for numeric diagonalization to calculate large 
number of energy eigenstates and eigenfunctions. 

Even when one uses such a sophisticated tool, the high- 
est eigenvalues cannot be used due to accumulation of 
numerical errors. In order to assess the quality of the 
obtained results for higher energy eigenstates, it is nec- 
essary to compare the numerical results with some known 
properties of the physical system. One such property is 



Psc(E) 



■ M 
,2irh 2 



r(f) 



d d x e(E-V(x)) (£7- V(x))* 



where O is the Heaviside step-function. For the quartic 
anharmonic potential (JSj) in d = 1 the density of states 
can be expressed in terms of the complete elliptic integral 
of the first kind K(k) = F(tt/2,A;) H, 



Psc(E) 



v /2M/ 7 r 2 fi 2 



(M 2 w 4 /4 + gE/Q) 1 / 4 




Mw 2 /4 



y/M 2 u A /A + gE/& 



(13) 



In practical applications, especially in d = 1, it might 
be difficult to compare directly semiclassical approxima- 
tion for density of states and numerically obtained his- 
togram for p(E), since energy levels are usually not de- 
generated, so the spectrum is very sparse. In order to 
have sufficient statistics for a reasonable histogram, one 
has to use large value for bin size, and effectively the 
whole numerically available spectrum is reduced to just 
a few bins. For this reason, it is more instructive to study 
the cumulative density of states, 



i(E) = f 

JVr, 



dE' p(E') , 



(14) 



which counts the number of energy eigenstates smaller 
or equal to E. For quartic anharmonic oscillator the cu- 
mulative density of states is given by the above integral 
of the complete elliptic integral of the first kind, and can 
be calculated numerically. Fig. [3] gives comparison of cu- 
mulative density of states calculated from our numerical 
diagonalization results and semiclassical approximation 
n sc (E). As expected, the agreement is excellent up to 



6 



300 
250 
200 
§ 150 

K 

100 
50 




g= 12,M = -10,L/A = 8192 
g = 12,M = -1,L/A = 8192 
g = 48,M = 1,L/A = 8192 
g = 48,M= 1,L/A= 100 



200 400 600 800 1000 1200 1400 1600 1800 



FIG. 3: (Color online) Cumulative distribution of the density 
of numerically obtained energy eigenstates for the quartic an- 
harmonic and double-well potential, for w = 1, t = 0.02, 
p = 21 and the following values of diagonalization parame- 
ters: L — 10 for <; = 12 and L = 8 for g = 48. The discretiza- 
tion step is given on the graph by the value of L/A, top to 
bottom. Long-dashed lines give corresponding semiclassical 
approximations for the cumulative density of states. 



very high values of energies, where numerical diagonal- 
ization eventually fails due to the finite number of cal- 
culated energy eigenvalues and effects of discretization. 
Such behavior can be improved by using finer discretiza- 
tion (smaller spacing size) , as illustrated by two different 
discretization steps for g = 48, M = 1 in Fig. [3J Such 
analysis can be used to assess the obtained spectrum 
and determine the number of reliable energy eigenval- 
ues. Typically we can achieve up to 10 4 reliable energy 
eigenlevels with simulations on a single CPU. Note that 
the computer double precision accuracy of 10 -16 imposes 
the limit on the maximal accessible energy eigenvalue 
e ~tE max „ 10 -i 6) i e Emax „ (i6i gio)/i. For Fig. E] 

we get E max ~ 1840, which is above the limit imposed 
by the discretization used, as we can see from the graph. 

In order to further demonstrate the applicability of 
the method, we also present numerical results for the 
modified Poschl- Teller model 

a 2 A(A-l) 
2 



V(x) 



(15) 



cosh ax 

which has only a finite set of discrete energy eigenlevels 
Ek = — a 2 (A — 1 — k) 2 /2 for integer k from the interval 
< k < A — 1. Energy eigenvalues and eigenfunctions 
of this model are analytically known, and we will use 
them to further test our method. Effective actions to 
very high order are available also for this potential [2l| . 
and we use them for numerical diagonalization of the 
evolution operator. Naturally, the diagonalization will 
give as many eigenvalues and eigenvectors as the size of 
the matrix S, but only the first few can be interpreted 
as bound states of the potential, according to the above 
condition < k < A — 1. 
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FIG. 4: ((Color online) Deviations |^ p) (A,L,t) - E% xact \ as 
a function of L for (a) k = and (b) k = 5, and (c) as a 
function of t for k = 0, for the modified Poschl- Teller poten- 
tial. Energy eigenvalues are obtained using effective action 
levels p = 1, 5, 9, 13, 17, 21 and t — 0.1 in graphs (a) and (b), 
and p = 1,3, 5, 7, 9, 11, 13 and L — 5 in graph (c), with the 
parameters a = 0.5, A = 15.5, A = 0.02. Dashed lines in the 
graph (c) correspond to the known discretization error 



Fig. gives the analysis of errors in the ground en- 
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ergy due to the space cutoff, while Fig. [4)3 gives the corre- 
sponding analysis of L-errors for numerical calculation of 
the energy level E$. As we can see, the behavior of errors 
is the same as for the case of anharmonic oscillator, and 
we are again able to obtain high accuracy results. Fig.|4j; 
gives the time dependence of errors in ground energy ob- 
tained by numerical diagonalization using different levels 
p of effective actions. The scaling of errors proportional 
to t p is evident from the graph, as well as the discretiza- 
tion errors due to the finite discretization step A. To 
ensure that the effective potential is bounded from be- 
low, in this case we have to remove higher-order powers 
of discretized velocity 6 from the effective potential near 
x = 0, since such terms have non- vanishing negative coef- 
ficients in the vicinity of x = 0, due to a peculiar nature 
of the potential. In practical applications, one can use 
e.g. p = 1 effective action (which does not depend on 
8) near x = 0. As can be seen, this does not affect the 
obtained numerical results. 

Table IHIa gives the obtained energy spectra for the 
modified Poschl- Teller potential with the parameters a — 
0.5, A = 15.5. If necessary, the precision of obtained 



energy levels can be further increased by appropriately 
changing the discretization parameters. Contrary to the 
situation for anharmonic oscillator, where relative error 
of numerically calculated low-lying energy levels did not 
change significantly, here we see that the increase in the 
error is substantial. This is caused by the fact that this 
potential has only a small finite set of discrete bound 
states, so energy levels k ~ 10 correspond to the very top 
of the discrete spectrum. In practical applications such 
pathological situations are not encountered, but as we 
can see, even this can be dealt with by the proper choice 
of discretization parameters. The quality of numerically 
calculated eigenf unctions is assessed in Table IPTIb , where 
we give a symmetric matrix of scalar products (ipk \tpf xact ) 
of numerically calculated and analytic eigenfunctions. As 
we can see, the overlap between analytic and numeric 
eigenfunctions is excellent, and they are orthogonal with 
high precision, which is preserved even for higher energy 
levels. We have also verified that for parameters given in 
the caption of Table IIIII and with the discretization step 
of the order A = 10~ 3 eigenfunctions of all bound states 
can be accurately reproduced. 
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TABLE III: (a) Low-lying energy levels of the modified Poschl- Teller potential, obtained by diagonalization using level p = 21 
effective action. The parameters used: a = 0.5, A = 15.5, L — 5, A = 0.02, t = 0.1. (b) Symmetric table of scalar products 
{V'klV'f" ) of numerically calculated and analytic eigenstates for k, I = 0, 1, 2, 3, 4, 5. 



IV. NUMERICAL RESULTS FOR d = 2 MODELS 

In this section we illustrate the application of the nu- 
merical method based on the diagonalization of transition 
amplitudes on two d — 2 models. The first model is the 
anharmonic oscillator 

V(x, y) = \m ( W i - tf ) (x* + jf) + ±(x* + y 2 f 
= IjtfodCl - r 2 ) (x 2 + y 2 ) + A(* 2 + y 2 ) 2 , (16) 



which is used for a description of the trapping poten- 
tial used in a recent experiment with fast-rota ting Bosc- 
Einstein condensate of 87 Rb atoms [H, H3, The 
harmonic frequency of the trapping potential is partially 
compensated by the rotation frequency f2. The small 
quartic anharmonicity is used in order to allow the con- 
densate to be rotated at the critical frequency f2 = uj±, 
and still to remain confined. The ratio r = Q/uj± is 
used to express rotation frequency in suitable units of 
harmonic frequency uj±. 

The typical values of parameters used in the experi- 



8 




FIG. 5: (Color online) (a) Deviations from the ground energy \Eq (A,L,t) — Eo\ as a function of the space cutoff L and (b) 
as a function of the time t for a critically rotating gas of 87 Rb atoms in a d = 2 anharmonic trap with g ~ 2 ■ W [i g exp . The 
discretization parameters are A = 0.2, t — 0.1 on the graph (a), and L — 3.2 on the graph (b). Deviations are calculated using 
the ground energy Eq — 1.47714975357799(4) obtained with p — 21 effective action. The dashed line in graph (b) corresponds 
to the known discretization error 
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-1.0311383813261 


4.6528451852013 
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4.0097032385903 


-1.0311383813261 


6.2704552903671 
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4.0097032385903 


-0.95910186300510 


6.2704552903671 
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4.0116368851078 


-0.95910186300478 


6.7589882491411 
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4.0116368851078 


-0.86968170695135 


6.7589882491412 



TABLE IV: Low-lying energy levels of a rotating gas of Rb atoms in a d — 2 anharmonic trap, obtained using the level p — 21 
effective action. The discretization parameters are L = 14, A = 0.14, and t — 0.2. 



merit are wj_ = 2ir ■ 64.8 Hz and g — g exp — 1.56 • 10~ 10 
J/m . In Fig. [5] we have used much larger quartic cou- 
pling g w 2-10 3 g exp in order to increase the non-harmonic 
effects of the potential. Also, the graphs in Fig. are 
calculated for critical rotation (r = 1), where the poten- 
tial is reduced to a pure quartic interaction. The anal- 
ysis of errors is very similar as in the one-dimensional 
cases we studied in the previous section. The depen- 
dence of ground energy errors on the space cutoff L is 
shown in Fig. [5^,, and we see the usual saturation of er- 
rors for sufficiently large values of L. The saturated value 
rapidly decreases (by several orders of magnitude) as we 
increase the level p of the effective action used to calcu- 
late space-discretized matrix of the evolution operator. 
Fig. [SJd shows the time dependence of ground energy er- 
rors, which are found to fully agree with the scaling law t p 
for sufficiently fine discretization. Again, the discretiza- 
tion errors basically conform to the universal dependence 
given in Eq. (17) of our previous paper [l[. The dimen- 
sionality of the system introduces an overall multiplica- 



tive factor of 2, and the additional factor of 2 in the cosh 
term. 

Table HVl gives the numerically obtained energy eigen- 
values for different sets of parameters of the potential: 
non-rotating system, system with overcritical rotation 
(r = 1.05), and system with overcritical rotation, but 
with significantly larger anharmonicity (g — 10 3 g exp ). 
From the analysis of discretization errors and errors re- 
lated to the use of a chosen effective action level p, we can 
estimate the errors in found energy eigenvalues to be of 
the order 10 -15 , where we express energy in units of hu>±. 
The results in the Table HVl are obtained by numerical di- 
agonalization based on the C SPEEDUP code [2l| and 
the use of the LAPACK 26] library. The estimated error 
in energy eigenvalues is smaller than the (relative) error 
which can be achieved in typical C simulations, which is 
of the order 1CP 14 . This is easily verified, since for several 
different values of discretization parameters we get the 
same stable results shown in the table. Therefore, this 
table gives certain digits in all energy eigenvalues, and 
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the error can be cited as implicit (half of the last digit). 
This is good example for practical applications, where we 
have managed to eliminate all types of errors below the 
limit that can be seen due to inherent numerical errors 
of computer simulation. However, if such complete elim- 
ination of errors is not possible due to the limitations in 
computer memory or computation time, the analysis of 
errors presented in Fig. O allows us to reliably estimate 
numerical errors in energy eigenvalues. 

Fig. [5] shows the numerically obtained ground state for 
this two-dimensional potential for the case of overcriti- 
cal rotation. The ground state has the expected Mex- 
ican hat shape. The figure gives a three-dimensional 
plot of the ground state on the left, and the correspond- 
ing density plot on the right, with values of the wave 
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function mapped to colors. Fig. gives density plots of 
k = 1, 2, 3, 4 eigenfunctions for the same values of param- 
eters. The discretization is sufficiently fine (A = 0.25) in 
rescaled dimensionless units) so that all features of cal- 
culated eigenfunctions are clearly visible. 

The numerical study of this example related to Bose- 
Einstein condensation is chosen as an example where 
ground state eigenfunction is necessary with high reso- 
lution in order to calculate e.g. time-of-flight absorption 
graphs [32| and to study formation and evolution of vor- 
tices in the condensate. In addition to this, large num- 
bers of accurate eigenstates are needed for calculation of 
the condensation temperature, condensate fraction, and 
other static and dynamic properties of the condensate. 
For this reason, it is necessary to assess numerically ob- 
tained eigenstates and use only reliable ones in further 
calculation. As in the one-dimensional case, we will cal- 
culate the density of states p sc (E) in semiclassical ap- 
proximation, and use it as a criterion for the reliability 
of high-energy eigenstates. In d = 2, the density of states 
is given by a simple formula 



Pi 



(17) 



For the quartic anharmonic potential (|16p the density of 
states can be analytically calculated 



Psc(E) 




QMoj\ (1 



6M4 (1 - r 2 ) 



2AE 
9 



(18) 



FIG. 7: (Color online) Density plots of level k = 1,2,3,4 
eigenstates of a rotating gas of 87 Rb atoms in a d — 2 an- 
harmonic trap obtained using p — 21 effective action. The 
parameters are r = 1.05, g = g exp , L — 20, A = 0.25, t = 0.2. 



or, in dimensionless units used in all numerical calcula- 
tions, 
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FIG. 8: (Color online) (a) Distribution of the density of nu- 
merically obtained energy eigenstates and (b) cumulative dis- 
tribution of the density of numerically obtained energy eigen- 
states for non-rotating gas of 87 Rb atoms in a d = 2 anhar- 
monic trap, calculated with the level p — 21 effective action. 
The parameters are r = 0, g = g eX p, t = 0.2, while discretiza- 
tion parameters are given on the graph, corresponding to the 
curves top to bottom. Long-dashed lines on both graphs give 
the corresponding semiclassical approximations. 
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TABLE V: Parameters of the sextic potential ([20}. 

from Fig. [5^, and again sets the same limit of reliable 
energy levels for chosen discretization parameters. 

The second two-dimensional model we have studied 
numerically is a sextic anharmonic oscillator, 



V{x,y) = V x (x) + V y {y) + V xy (x - y) . 



(20) 



where Vi(x) = Vm (aiX 2 + bix 4 + ax 6 ). The values of 
the coefficients used are given in Table [Vl The study 
of this potential is motivated by Ref. |33|, where it has 
been used to investigate the transition from regular to 
chaotic classical motion. Fig. [5] shows the numerically 
obtained ground state for this two-dimensional potential, 
as a three-dimensional plot on the left, and as a density 
plot on the right. Fig. [TU] gives density plots of k = 
1, 3, 7, 8 eigenfunctions for the same values of parameters. 
The discretization is sufficiently fine (A = 0.04) so that 
we can resolve all details in the presented eigenstates. 
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Fig. [8^ shows the comparison of semiclassical approx- 
imation for the density of states, and the histogram for 
numerically obtained energy eigenvalues of the potential 
(p~6|) . Due to the high degeneracy of energy eigenstates in 
d = 2, the histogram of numerically found energy levels 
contains enough statistics over the whole region of en- 
ergies, and therefore can be used for assessment of the 
quality of numerical spectra. As we see, the agreement is 
better and better when we use finer space discretization. 
Depending on the needed number of energy levels and 
maximal value of the energy considered to be relevant for 
the calculation we can choose appropriate values of dis- 
cretization parameters that will provide reliable numeri- 
cal results up to desired energy value. For example, for 
the choice of discretization parameters L = 14, A = 0.14, 
we can reliably use energy levels up to E w 120 ftc^. 

Fig. [5Jd shows the comparison of cumulative density 
of states n(E) calculated for numerically obtained re- 
sults and in semiclassical approximation, by integrating 
the expression ([TO]) , which can be calculated analytically. 
The comparison of numerical and semiclassical cumula- 
tive density of states in Fig. [SJd verifies our conclusions 
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FIG. 9: (Color online) Ground state (as a 3-D and as a density 
plot) of a sextic anharmonic potential, obtained by diagonal- 
ization using the level p = 21 effective action. The parameters 
of the potential are given in the text. The diagonalization pa- 
rameters: L = 4, A = 0.04, t = 0.01. 
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FIG. 10: (Color online) Density plots of level k = 1,3,7,8 
eigenstates of a sextic anharmonic potential, obtained by di- 
agonalization using the level p = 21 effective action. The 
parameters of the potential are given in the text. The diago- 
nalization parameters: L — 4, A = 0.04, t = 0.01. 

We have demonstrated that the presented approach 
can be successfully used for numerical studies of lower- 
dimensional models. Note that in d = 3 the complexity 
of the algorithm and sizes of matrices to be diagonal- 
ized may practically limit the applicability to the cal- 
culation of only low- lying energy levels. Also, in this 
case it might be difficult to numerically obtain three- 
dimensional eigenfunctions on finer grids, since even 
moderate grids with 50-100 points in one dimension 
would require exact diagonalization of extremely large 
matrices. 

At the end, let us compare the complexity of the pre- 
sented approach and direct diagonalization of the space- 
discretized Hamiltonian, as well as finite-element meth- 
ods. The main difference in the complexity of algorithms 
is related to the exponential growth in the size of analytic 
expressions for the effective potential with the increase of 
the level p, as discussed in Ref. jl3j |. Therefore, the re- 
quired CPU time for construction of the matrix to be di- 
agonalized in the presented approach grows exponentially 
with the level p, while in other methods the construction 
of such a matrix does not require a significant amount 
of time. However, the time for exact diagonalization far 
outweighs the time needed for construction of even large 
matrices with moderate levels p of the order 10-20. The 
significant benefit of practically eliminating errors associ- 
ated with the time of propagation therefore fully justifies 
the use of the effective action approach. Of course, in 
practical applications one has to study the complexity of 
the algorithm and to choose the optimal level p which 
will sufficiently reduce the errors, while keeping the com- 
plexity of the calculation on the acceptable level. 



V. CONCLUSIONS 

In this paper we have presented a substantial improve- 
ment of previously introduced method Q for study or 
properties of quantum systems using numeric diagonal- 
ization of the space-discretized evolution operator. This 
approach allows exact numeric calculation of a large num- 
ber of energy eigenvalues and eigenstates of the system. 
Our previous paper [l| has presented detailed analysis of 
all types of discretization errors inherent to this method, 
which were not analyzed completely before. 

This paper resolves a key problem in practical appli- 
cations of this approach: accurate calculation of transi- 
tion amplitudes, matrix elements of the space-discretized 
evolution operator. Using recently introduced effective 
action approach [l3| that gives systematic short-time ex- 
pansion of the evolution operator, we can analytically 
calculate matrix elements of the evolution operator with 
high precision. This enables high precision calculation of 
energy eigenvalues and eigenstates, as was shown in this 
paper. 

The derived analytical estimates for all types of errors, 
including errors due to the approximative calculation of 
transition amplitudes, provide us with a way to choose 
optimal discretization parameters and to reduce overall 
errors in energy eigenvalues and eigenstates for many or- 
ders of magnitude, as was demonstrated for several one- 
and two-dimensional models. We have shown that nu- 
merical diagonalization of the space-discretized evolution 
operator can be successfully applied for studies of many 
interesting lower dimensional models. Due to the su- 
perior behavior of discretization and other errors in this 
method compared to methods based on diagonalization of 
the discretized Hamilton operator and related methods, 
the presented approach is a method of choice for numer- 
ical studies of lower-dimensional physical systems. The 
authors are already using this approach for numerical in- 
vestigation of properties of fast rotating Bose-Einstein 
condensates [32j . and plan to use it for the treatment of 
dilute quantum gases in a disordered environment. An- 
other interesting line of research would be combining the 
present method with the density matrix renormalization 
group (DMRG) approach @,[34|. 



Acknowledgments 

This work was supported in part by the Ministry of Sci- 
ence and Technological Development of the Republic of 
Serbia, under project No. OI141035 and bilateral project 
PI-BEC funded jointly with the German Academic Ex- 
change Service (DAAD), and the European Commission 
under EU Centre of Excellence grant CX-CMCS. Numer- 
ical simulations were run on the AEGIS e-Infrastructure, 
supported in part by FP7 projects EGEE-III and SEE- 
GRID-SCI. 



12 



I. Vidanovic, A. Bogojevic, A. Belie, Phys. Rev. E 80 [18 
(2009) 066705; larXiv:0911.5145l [19 
A. Sethia, S. Sanyal, Y. Singh, J. Chem. Phys. 93 (1990) 
7268. [20 
A. Sethia, S. Sanyal, F. Hirata, Chem. Phys. Lett. 315 
(1999) 299. [21 
A. Sethia, S. Sanyal, F. Hirata, J. Chem. Phys. 114 
(2001) 5097. [22 
S. Sanyal, A. Sethia, Chem. Phys. Lett. 404 (2005) 192. [23 
T. L. Beck, Rev. Mod. Phys. 72 (2000) 1041. 
M. A. Martin-Delgado, G. Sierra, R. M. Noack, J. Phys. [24 
A 32 (1999) 6079. [25 
J. R. Chelikowsky, N. Troullier, K. Wu, Y. Saad, Phys. 
Rev. B 50 (1994) 11355. [26 
P. Maragakis, J. M. Soler, E. Kaxiras, Phys. Rev. B 64 
(2001) 193101. [27 
S. A. Chin, S. Janecek, E. Krotscheck, Comp. Phys. 
Comm. 180 (2009) 1700. 
S. Janecek, E. Krotscheck, Comp. Phys. Comm. 178 [28 
(2008) 835. 

M. Aichinger, S. A. Chin, E. Krotscheck, Comp. Phys. 
Comm. 171 (2005) 197. [29 
A. Balaz, A. Bogojevic, I. Vidanovic, A. Pelster, Phys. 
Rev. E 79 (2009) 036701. [30 
A. Bogojevic, A. Balaz, and A. Belie, Phys. Rev. Lett. 
94 (2005) 180403. [31 
A. Bogojevic, A. Balaz, and A. Belie, Phys. Rev. B 72 [32 
(2005) 064302. 

A. Bogojevic, A. Balaz, and A. Belie, Phys. Lett. A 344 [33 
(2005) 84. 

A. Bogojevic, I. Vidanovic, A. Balaz, A. Belie, Phys. [34 
Lett. A 372 (2008) 3341. 



D. M. Ceperley, Rev. Mod. Phys. 67 (1995) 279. 

M. Boninsegni, N. Prokof'ev, B. Svistunov, Phys. Rev. 

Lett. 96 (2006) 070601. 

A. Bogojevic, A. Balaz, A. Belie, Phys. Rev. E 72 (2005) 
036128. 

SPEEDUP C and Mathematica code, 
http : //www . scl . rs/speedup/ 

W. Janke, H. Kleinert, Phys. Rev. Lett. 75 (1995) 2787. 

G. Barton, A. J. Bray, A. J. McKane, Am. J. Phys. 58 
(1990) 751. 

D. H. Berman, Am. J. Phys. 59 (1991) 937. 
Mathematica symbolic calculation software package, 
http : //www . wolfr am. com/| 

LAPACK Linear A lgebra PACKage, 

http : / /www . netlib . org/lapack/ 

H. Kleinert, Path Integrals in Quantum Mechanics, 
Statistics, Polymer Physics, and Financial Markets, 4th 
edition, World Scientific, Singapore, 2006. 

I. S. Gradshteyn, I. M. Ryzhik, Table of integrals, series, 
and products, 6th edition, Academic Press, San Diego, 
2000. 

V. Bretin, S. Stock, Y. Seurin, J. Dalibard, Phys. Rev. 
Lett. 92 (2004) 050403. 

S. Stock, B. Battelier, V. Bretin, Z. Hadzibabic, J. Dal- 
ibard, Laser Phys. Lett. 2 (2005) 275. 
S. Kling, A. Pelster, Phys. Rev. A 76 (2007) 023609. 
A. Balaz, I. Vidanovic, A. Bogojevic, A. Pelster, in prepa- 
ration. 

T. H. Seligman, J. J. M. Verbaarschot, M. R. Zirnbauer, 

Phys. Rev. Lett. 53 (1984) 215. 

U. Schollwock, Rev. Mod. Phys. 77 (2005) 259. 



